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ABSTRACT 

We use single-dish radio spectra of known 22 GHz H 2 O megamasers, primarily gathered from 
the large dataset observed by the Megamaser Cosmology Project, to identify Keplerian accretion 
disks and to investigate several aspects of the disk physics. We test a mechanism for maser 
excitation proposed by Maoz & McKee (1998), whereby population inversion arises in gas behind 
spiral shocks traveling through the disk. Though the flux of redshifted features is larger on 
average than that of blueshifted features, in support of the model, the high-velocity features 
show none of the predicted systematic velocity drifts. We find rapid intra-day variability in the 
maser spectrum of ESO 558-G009 that is likely the result of interstellar scintillation, for which we 
favor a nearby (D ss 70 pc) scattering screen. In a search for reverberation in six well-sampled 
sources, we find that any radially-propagating signal must be contributing <10% of the total 
variability. We also set limits on the magnetic field strengths in seven sources, using strong 
flaring events to check for the presence of Zeeman splitting. These limits are typically 200-300 
mG (lcr), but our most stringent limits reach down to 73 rnG for the galaxy NGC 1194. 


Subject headings: accretion, accretion disks — magnetic fields — masers — galaxies: active — galaxies: 
nuclei — galaxies: Seyfert — radio lines: galaxies — scattering 


1. Introduction 

The Megamaser Cosmology Project (MCP) 
aims to determine the value of Hq by measur- 
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ing angular-diameter distances to galaxies in the 
Hubble flow. Using the megamaser technique pi¬ 
oneered on the galaxy NGC 4258 (see Herrnstein 
et al. 1999), the MCP has published distances to 
the galaxies UGC 3789 (Reid et al. 2013), NGC 
6264 (Kuo et al. 2013), and NGC 6323 (Kuo et al. 
2015), and additional galaxies are currently be¬ 
ing measured. The ongoing project is a multi¬ 
year effort of surveying, monitoring, and mapping 
maser disks using the Robert C. Byrd Green Bank 
Telescope (GBT), the Karl G. Jansky Very Large 
Array (VLA), and the Very Long Baseline Array 
(VLBA) plus the 100-meter Effelsberg telescope. 

The MCP’s monitoring campaign uses the GBT 
of the National Radio Astronomy Observatory 
(NRAO) to take regular (^monthly) spectra of 
megamaser sources targeted for distance measure¬ 
ments. These spectra are used to measure the 
accelerations of maser features as part of the de¬ 
termination of Hq. Here we take advantage of 
this rich dataset to probe the innermost parsec 
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(~0.1-0.5 pc) of the AGN. These size scales are 
roughly an order of magnitude smaller than the 
dust structures that have been resolved by opti¬ 
cal/infrared interferometric studies of the torus re¬ 
gion in nearby AGN (see, e.g., Jaffe et al. 2004). 

The structure of this paper is as follows. The 
observations and data reduction procedures are 
described in §2. We present the data in §3 and 
§5.1, in the form of time-averaged and dynamic 
spectra, respectively. In §4 we examine a theory 
of disk maser excitation proposed by Maoz & Mc¬ 
Kee (1998) (hereafter MM98), in §5.2 we present 
evidence for the presence of interstellar scintilla¬ 
tion in ESO 558-G009, in §6 we check the maser 
disks for signs of propagating disturbances, and in 
§7 we use the spectra to place limits on the mag¬ 
netic field strengths in the maser disks. 

2. Observations and data reduction 

The analyses presented in this paper are based 
on 22 GHz water maser spectra, almost all of 
which were taken using the GBT over the pe¬ 
riod March 2003 April 2015. The majority of 
these spectra were obtained as part of the survey 
and monitoring components of the MCP; see Reid 
et al. (2009) and Braatz et al. (2010) for details. 
We include several non-MCP spectra from the 
NRAO data archive, most notably for the galaxies 
NGC 4258 and NGC 3393. 

For each MCP spectrum the GBT spectrom¬ 
eter was configured with two 200 MHz spectral 
windows, one of which was centered on the reces¬ 
sion velocity of the galaxy while the other was off¬ 
set redward by 180 MHz. Each window had 8192 
channels spaced at 24 kHz channel width, which 
at 22 GHz corresponds to approximately 0.33 km 
s -1 . Both left circular polarization (LCP) and 
right circular polarization (RCP) were observed 
simultaneously in each of the two beams of the 
K-band receiver, and the telescope was nodded 
on a 2.5-minute cycle to alternate which beam 
was pointed at the target. Observations after 
May 2011 used two of the seven beams of the K- 
band Focal Plane Array (KFPA) in the same nod¬ 
ding scheme. Integration times for the monitored 
sources were typically between 1 and 3 hours dur¬ 
ing a single observing session. 

We reduced GBT data using the same methods 
outlined in previous MCP papers (see, e.g., Braatz 


et al. 2010). Our measurements of Zeeman split¬ 
ting (§7.2) use spectra at their native resolution, 
prior to Hanning smoothing. 

Integrated line fluxes in some of our spectra are 
affected by a broad (~1500 km s -1 ) sinusoidal 
baseline ripple. The baseline ripples are gener¬ 
ally comparable in amplitude to the RMS channel 
noise, but their contributions to the flux measure¬ 
ments can be the dominant source of uncertainty 
for our best-sampled sources. To characterize the 
flux uncertainty from the baseline ripple, we av¬ 
eraged the frequency-offset spectral windows from 
each observation. These spectra are free of maser 
emission, were taken concurrently with the science 
spectra using the same instrument configuration 
on the GBT, and have all undergone the same data 
reduction procedure. We measured the RMS of 
the integrated flux inside a boxcar window placed 
randomly inside the averaged spectrum, as a func¬ 
tion of the spectral width of that window. The line 
flux uncertainty behaves approximately quadrat- 
ically as a function of window width, reaching a 
maximum of ~0.1 Jy km s -1 for a window width 
of ~750 km s -1 . We thus assign a baseline ripple 
uncertainty to each line flux measurement that fol¬ 
lows the empirical relation given by 

0 ’S,o.i = — (Au ) 750 + 2 (Au ) 750 (1) 

Here, crs,o.i is the baseline rippler’s contribution 
to the line flux uncertainty (in units of 0.1 Jy km 
s^ 1 ) and (An )750 is the spectral window width (in 
units of 750 km s -1 ). 

3. Identifying Keplerian disk megamasers 

Our aim is to examine the spectral characteris¬ 
tics of maser emission from accretion disks, includ¬ 
ing flux ratios, secular velocity drifts, and variabil¬ 
ity. We thus seek to identify maser systems with 
spectra dominated by emission from edge-on, Ke¬ 
plerian disks. 

To date, 16 megamaser disk systems have pub¬ 
lished VLBI maps. Eight of these were mapped by 
the MCP (NGC 1194, NGC 2273, Mrk 1419, NGC 
4388, NGC 6323 in Kuo et al. 2011; UGC 3789 in 
Reid et al. 2009; NGC 6264 in Kuo et al. 2013; 
and NGC 5765b in Gao et al. submitted), and 
eight were mapped by other groups (NGC 1068 in 
Greenhill & Gwinn 1997; NGC 4945 in Greenhill 
et al. 1997b; NGC 5793 in Hagiwara et al. 2001; 
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Circinus in Greenhill et al. 2003b; NGC 3079 in 
Kondratko et al. 2005; NGC 4258 in Argon et al. 
2007; NGC 3393 in Kondratko et al. 2008; and IC 
1481 in Mamyoda et al. 2009). Of these 16 mapped 
disks, nine have “clean” Keplerian rotation curves, 
and all nine share a distinctive single-dish spectral 
profile. To maximize the uniformity and size of the 
sample for the analysis in this paper, we therefore 
selected sources based on the appearance of their 
single-dish (usually GBT) spectra. 

A “clean” disk nregamaser is an edge-on maser 
in Keplerian rotation around the central SMBH, 
in which the disk maser emission dominates over 
any jet or outflow maser components. These sys¬ 
tems have characteristic spectra that are marked 
by three sets of maser components. The “sys¬ 
temic” set of features coincides roughly with the 
recession velocity of the galaxy, and the masing 
arises along a line of sight through the disk to the 
central AGN. The two “high-velocity” sets of fea¬ 
tures (the “redshifted features” and “blueshifted 
features”) are spectrally offset to either side of the 
galaxy’s recession velocity. These features arise 
from the midline of the accretion disk, along lines 
of sight that are tangent to the orbital motion 
(which ensures velocity coherence throughout the 
column of gas). For an edge-on disk, the midline 
is the diameter through the disk that falls perpen¬ 
dicular to the line of sight. 

To select clean disk nregamasers, we use the 
following criteria. The spectra must show at least 
two of the three expected distinct sets of maser fea¬ 
tures (in maser disks with only two sets of features, 
the third set is presumably present but below the 
detection threshold). Furthermore, at least one of 
the sets of features should have components that 
are offset from the recession velocity by at least 
300 km s —1 (an empirically-determined high galac¬ 
tic rotation cutoff; see Cresci et al. 2009), to avoid 
contaminating the sample with interstellar masers 
(from, e.g., a strong starburst) and sub-Keplerian 
rotators. For a spectrum with only two sets of 
maser features, we require either that one of these 
feature sets be coincident with the recession veloc¬ 
ity of the galaxy or that both feature sets be offset 
from the recession velocity by at least 300 km s _1 . 

Though we have attempted to be comprehen¬ 
sive in our selection of sources, there are several 
known disk (or disk-like) H 2 O megamasers that 
do not make it into our sample because the disk 


emission is contaminated by non-disk components. 
Circinus (Gardner & Whiteoak 1982) contains a 
masing accretion disk, but it also has maser emis¬ 
sion associated with an outflow (Greenhill et al. 
2003b). Similarly, NGC 1068 (Claussen et al. 
1984) has maser emission arising from both a disk 
and a radio jet encountering a dense molecular 
cloud (Gallimore et al. 1996). Complexities like 
these confuse the maser spectrum and make it dif¬ 
ficult to associate individual spectral features with 
either the disk or outflow/jet components without 
a VLBI map. For this reason, none of these sources 
passes our selection criteria. 

The final list of 32 clean meganraser disk sys¬ 
tems used in our study is given in Table 1. Fig¬ 
ure 1 shows the weighted average spectra of these 
sources, where the weighting r/T S y S (where r is 
the exposure time and T sys is the system temper¬ 
ature) was chosen to minimize the RMS noise of 
each spectrum. The emission from the remaining 
~130 known water megamaser galaxies may arise 
from nuclear sources other than the accretion disk 
(e.g., molecular gas in an outflow) or from extranu- 
clear sources elsewhere within these galaxies (e.g., 
star-forming regions). 

For completeness, we reproduce the spectrum 
of ESO 269-G012 in Figure 1 from Greenhill et al. 
(2003a); see that paper for details about the ob¬ 
servation and data reduction. 

3.1. Observed properties of disk mega¬ 
masers 

Table 1 also lists several observational prop¬ 
erties of each galaxy. We obtained the reces¬ 
sion velocities from the NASA/IPAC Extragalac- 
tic Database (NED), favoring velocities measured 
from neutral hydrogen (HI) over those made from 
optical lines. HI measurements have the advan¬ 
tage that they average over all internal motions 
of a galaxy, while optical lines are preferentially 
emitted from regions with a sufficiently energetic 
radiation field to excite the transitions. In the case 
of active galaxies, like those present in our sample, 
the optical emission could very well be dominated 
by gas that is kinematically driven by the nuclear 
activity (e.g., outflows). This could result in a 
systematic offset between the recession velocity of 
the galaxy and the velocity measured using optical 
lines (e.g., Comerford et al. 2013). We do see such 
offsets in several of the spectra shown in Figure 1. 
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Table 1 
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Fig. 1. - Spectra for disk megamasers used in our analysis of the MM98 model. Each spectrum is a weighted 
average (see §3) taken over all epochs; the date of the first epoch is located at the top right. Galaxy recession 
velocities and associated la errors (see Table 1) are overplotted in red. 
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Fig. 1. — (continued) 



2008 Nov 09 


ilLkiiMriililiim MaikliLy 


14000 


14500 15000 15500 








0.015 

0.010 

0.005 

0.000 

-0.005 


NGC4388 


2003 Mar 07 




2000 2200 2400 2600 2800 3000 

Velocity (km/s) 




U.1Z 

0.10 

ESO 269-G012 + 2002 

Jun 1 

0.06 

L 1 


- 

0.04 




0.02 

0.00 

0.02 

ft 

hr 




4500 5000 5500 

Velocity (km/s) 


6 


































































Fig. 1. — (continued) 
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Fig. 1. — (continued) 
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We measured line fluxes separately for each set 
of features: blueshifted, systemic, and redshifted. 
To maximize the signal-to-noise for those spectra 
with weak features, we integrated only over spec¬ 
tral windows that contained clear signal. In some 
cases this meant integrating over several distinct, 
narrow windows to obtain the total line flux for 
that set of features. For several spectra, the sys¬ 
temic set of features is absent; in these cases we 
list an upper limit on the line flux for the systemic 
features obtained by integrating over the spectral 
region located between the high-velocity features 
(i.e., the region redward of the blueshifted features 
and blueward of the redshifted features). 

To obtain the total isotropic luminosities listed 
in Table 1, we integrated each spectrum across the 
full span of maser emission. For a measured line 
flux S, the isotropic luminosity is given by 


Here, v is the recession velocity of the galaxy. This 
expression is accurate for low-redsliift (z < 0.1) 
sources, and all of our galaxies fall into this cat¬ 
egory so we use it throughout. In our calcula¬ 
tions, we assume a Hubble constant of Hq = 70 
km s -1 Mpc -1 . Figure 2 shows a histogram of 
the isotropic luminosities. To alleviate the some¬ 
what arbitrary nature of histogram bin sizes and 
endpoints, we have overplotted a kernel density es¬ 
timate using a Gaussian kernel. The area of each 
kernel is equal to that of a histogram bin with 
a bin size determined using Silverman’s rule; see 
Appendix A for details. 

The measured isotropic luminosities span over 
two orders of magnitude, and the observed dis¬ 
tribution (see Figure 2) appears to be consistent 
with a sensitivity-limited sample (i.e, the high¬ 
est luminosity masers tend to be found at large 
distances, and vice versa). While some of this 
spread is undoubtedly caused by intrinsic power 
differences among the many systems, most of it is 
likely the result of viewing angles. Though the ex¬ 
act angular dependence of the maser emission is a 
strong function of the source geometry and satu¬ 
ration, it always drops off exponentially from the 
beam center, which falls along the path of maxi¬ 
mum gain. (For an in-depth discussion of maser 
beaming, we refer the interested reader to Elitzur 
1992.) Thus, even a slight (<5°) inclination of the 


maser beam from the line of sight could cause the 
observed intensity to drop by an order of magni¬ 
tude or more. This is especially true if the masers 
are unsaturated. The unknown contribution from 
maser beaming precludes us from correcting the 
Malmquist bias and turning Figure 2 into a true 
luminosity function. 

4. Testing a model of disk maser excita¬ 
tion 

In their 1998 paper, Maoz & McKee (MM98) 
sought to explain the observation in NGC 4258 
that the line flux of the redshifted features is much 
higher than the line flux of the blueshifted fea¬ 
tures. In their model, population inversion (and 
thus masing) only occurs in post-shock gas on the 
trailing edge of a spiral shock in the accretion disk. 
Observed high-velocity maser features then occur 
wherever the line of sight falls tangent to a shock 
front, for an edge-on disk system. 

The geometry of the trailing spiral shocks 
causes redshifted maser emission to preferentially 
originate from the region of the disk that lies 
in front of the midline, while blueshifted maser 
emission arises from behind the midline. The 
blueshifted photons would thus pass through a 
sightline of velocity-coherent (but noninverted) 
gas, leading to absorption that is not present for 
the redshifted photons. The model thereby pre¬ 
dicts that the redshifted high-velocity features ob¬ 
served for disk maser systems should be systemat¬ 
ically stronger than the blueshifted high-velocity 
features. See Fig. 1 in MM98 for an illustration 
of this geometry. 

Owing to their offsets from the midline, the 
MM98 model predicts nonzero line-of-sight “ac¬ 
celerations”; specifically, the blueshifted features 
should show a mean positive acceleration while the 
redshifted features show a negative one. These 
arise because as the trailing spiral shock passes 
through the disk, the inversion region (and thus 
the segment of spiral structure that is tangent to 
the line of sight) moves radially outwards with 
time. The line-of-sight component of the veloc¬ 
ity decreases in magnitude with increasing radius, 
so the result is an observed velocity drift in the 
high-velocity maser lines. Though such behavior 
mimics an acceleration, it is actually tracing the 
rotating spiral structure rather than the Keple- 
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Fig. 2.— Histogram showing the distribution of isotropic luminosities for our sample of disk megamasers. 
The solid black line shows the kernel density estimation obtained using a normal kernel, with Silverman’s 
rule applied for the kernel Gaussian width (see Appendix A). 
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rian motion of the gas in the disk, and we there¬ 
fore refer to the phenomenon as a “velocity drift” 
rather than as an acceleration (see §4.2 for de¬ 
tails). This prediction runs counter to that of 
the “standard” model, which has an entirely mas- 
ing disk with high-velocity features falling close 
to the midline. The standard model predicts that 
the high-velocity features should have nearly zero 
line-of-sight accelerations on average. 

The model proposed by MM98 was inspired by 
the red-blue flux asymmetry in NGC 4258, which 
we note from Table 1 has a uniquely high value of 
log(-R) = +1.42 not seen in any other maser disk. 
It is an open question whether such an excitation 
mechanism applies to maser disks in general; in¬ 
deed, it is an open question whether this mech¬ 
anism even holds for NGC 4258 (see, e.g., Bragg 
et al. 2000). We checked this model by measuring 
the flux asymmetry and velocity drifts of high- 
velocity features in our Keplerian disk sample. 

4.1. Statistical analysis 

For each disk maser in our sample we made 
a weighted average spectrum from all epochs of 
observation (see Figure 1 and §3). The averag¬ 
ing reduces the noise and mitigates the effects 
of variability. We then identified the regions of 
each spectrum corresponding to the redshifted and 
blueshifted high-velocity features. By integrating 
over these spectral segments, we obtained the red- 
shifted and blueshifted fluxes. The ratio, R , of the 
redshifted to the blueshifted flux should be greater 
than 1 for the MM98 model. The values of log (R) 
for our sample are listed in Table 1 and their his¬ 
togram is plotted in Figure 3. 

The null hypothesis is that the redshifted and 
blueshifted fluxes are on average equal; that is, 
the logarithm of the ratio of the redshifted to the 
blueshifted flux should be a distribution centered 
on zero. We use the logarithm of the flux ratios 
(rather than the ratios themselves) to avoid the 
skewing of the distribution that arises from a di¬ 
rect ratio. 

To test whether our results are consistent with 
the null hypothesis, we employ a likelihood anal¬ 
ysis to determine whether the sample we observe 
has been drawn from a parent population with an 
intrinsic flux ratio distribution centered on zero. 
The data point corresponding to NGC 4258 is not 


included in this analysis, as it was used to gen¬ 
erate the original hypothesis. Here we utilize a 
technique analogous to that presented in Richards 
et al. (2011). 

To simplify notation, we define X = log (R), 
where R = p/P is the ratio of the redshifted flux 
(denoted p) to the blueshifted flux (denoted P). 
We assume that the parent distribution of X is 
a Gaussian centered on Xo, with a standard de¬ 
viation of (Jo. We also assume that the observa¬ 
tional uncertainties associated with each measure¬ 
ment are normally distributed about the intrinsic 
value for that measurement. 

For a single observation of a source with intrin¬ 
sic redshifted flux of p t , the probability to observe 
the value pi with uncertainty c+j is given by 


Pr = 


1 



0 pt - Pi) 2 

Zvr,i 


( 3 ) 


Similarly for an observation of a source with in¬ 
trinsic blueshifted flux of p t , the probability to 
observe the value Pi with uncertainty o+j will be 


P b = 


1 



(Pt - Pi ) 2 

2c7 M 


( 4 ) 


We also have the probability for the source to 
have an intrinsic flux ratio of X t = log(p t //3 t ), 
given the parent distribution 


Pt = 



(x t -x 0 f - 

2cr o 


( 5 ) 


The resulting likelihood of the observation is then 
given by an integral over the product of these 
probability density functions, 


t,. = I / P r P b Ptdptdp t . (6) 
Jo Jo 

For N observations, the joint likelihood will then 
be the product of the individual measurement like¬ 
lihoods: 


N 

C(X 0l a 0 ) (7) 

2—1 

Once the joint likelihood function is known, 
we can marginalize over the parameter a 0 . The 
marginalized likelihood, £(X 0 ) (shown in Figure 
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log(R) 


Fig. 3.— Histogram showing the distribution of the logarithm of the red/blue flux ratios for our sample of 
disk megamasers. The solid black line again shows the kernel density estimation, obtained using the same 
normalization as in Figure 2. NGC 4258 occupies the rightmost histogram bin, causing the red tail of the 
distribution to be noticeably longer and heavier than the blue tail. Though we include it in this plot, NGC 
4258 was not included in the statistical analysis performed in §4.1 to avoid biasing the results (i.e., since the 
proposed hypothesis was based on observations of NGC 4258, its observed properties necessarily agree with 
the hypothesis). 
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Fig. 4.— The normalized likelihood for the model 
presented in §4.1 as a function of the average flux 
ratio Xq = (log(i?)), marginalized over ao- Ranges 
corresponding to lcr and 2cr are shown as dashed 
and dotted vertical lines, respectively. 


4), can then be integrated to determine the frac¬ 
tion of the likelihood that falls below Xq = 0: 


P = 


C(X 0 )dX o 


C(X 0 )dX 0 


( 8 ) 

Evaluating p for the flux values listed in Table 1 
yields p = 0.020. The likelihood analysis therefore 
rejects the null hypothesis at the 2cr level. 


NGC 4258 stands out as an 18<r outlier, which 
most likely indicates that this Gaussian model is 
not a good description of the parent population. 
Nevertheless, it is sufficient to show that the null 
hypothesis is at least moderately discrepant with 
the data and that NGC 4258 is substantially re¬ 
moved from the bulk of the observed distribution. 


4.2. Velocity drifts of high-velocity fea¬ 
tures 

The MM98 model also predicts that the high- 
velocity maser features will be systematically off¬ 
set from the midline of the disk, and that they 


should thus exhibit nonzero line-of-sight velocity 
drifts as the spiral structure rotates. For spiral 
shocks having a pitch angle of 9 p (the pitch angle 
is the opening angle of the spiral, defined at any 
point to be the complement of the angle between 
the tangent to the spiral and the outward radial 
direction from the black hole), we can calculate a 
characteristic value for the velocity drifts expected 
for the high-velocity features. For a logarithmic 
spiral, the MM98 model predicts a velocity drift 
of 


1*1 = °-°5 (^y) km s 1 yr 1 . (9) 

This drift is towards smaller rotation velocities, 
and it is shared by all high-velocity masers. The 
observed velocity drift, in this model, is caused by 
the passage of the trailing spiral structure through 
the gas disk; it is not a centripetal acceleration 
from the Keplerian rotation of the gas. As the 
spiral shock moves through the disk, the portion 
tangent to the line of sight intercepts gas farther 
out in radius, which has a lower rotational velocity. 
Thus we would expect to observe a negative line- 
of-sight velocity drift for the redshifted features 
and a positive drift for the blueslrifted features. 

Bragg et al. (2000) measured velocity drifts in 
NGC 4258, and showed that the values were incon¬ 
sistent with the predictions of the MM98 model. 
They established that no choice of pitch angle can 
reproduce their data, as statistically significant 
measurements of both negative and positive ve¬ 
locity drifts were made for both sets of features. 
These results were corroborated by Humphreys 
et al. (2008), who used an increased number of 
epochs to further refine the measurements. Ta¬ 
ble 2 lists published measurements of high-velocity 
drifts for several other megamaser disks, plus our 
new measurements, where we have estimated ve¬ 
locity drifts for several additional galaxies using 
the eye-tracking method described in Kuo et al. 
(2013). To account for systematic uncertainties, 
we also adopt the error floor of 0.3 km s -1 yr -1 
from Kuo et al. (2013) for all new acceleration 
measurements. 

Nine of the 22 velocity drift measurements 
(counting redshifted and blueslrifted separately) 
presented in Table 2 are incompatible with the 
MM98 model (i.e., negative blueslrifted velocity 
drifts or positive redshifted velocity drifts). For 
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those values which are compatible, we used Equa¬ 
tion 9 to assign a maximum pitch angle of any 
spiral structure that is consistent with the mea¬ 
sured drifts. As a comparison, the minimum pitch 
angle in NGC 4258 (obtained by assuming that 
the spatial grouping of the blueshifted features 
arise from consecutive windings of a single log¬ 
arithmic spiral) is about 9 p > 1.7° (Humphreys 
et al. 2013). 

4.3. Discussion 

Our analysis of the flux ratio data indicate a 
small deviation from the null hypothesis, in favor 
of the MM98 model. However, the measured ve¬ 
locity drifts of high-velocity features do not match 
the MM98 predictions (Table 2). The maser fea¬ 
tures are equally likely to have a positive drift 
as a negative one, regardless of whether they’re 
blueshifted or redshifted (6 of 11 targets display 
negative velocity drifts for both sets of features). 
Furthermore, though we have reported only the 
averaged values for the redshifted and blueshifted 
velocity drifts for each target, several of these tar¬ 
gets have statistically significant measurements of 
both negative and positive drifts within the same 
set of features. On the whole, the high-velocity 
drifts are consistent with masing gas that is near 
the midline of the disk (i.e., any observed velocity 
drifts can be explained as centripetal accelerations 
caused by small offsets on both sides of the mid¬ 
line). 

We note that the MM98 model is based on 
the characteristics of NGC 4258, which has an 
atypically large flux ratio between the redshifted 
and blueshifted high-velocity features. This ap¬ 
parent anomaly could be the result of a selection 
bias. If NGC 4258 were located at a distance of 
~ 100 Mpc, which is more typical of our sample, 
it would likely not have been identified as a disk 
maser. The systemic features would peak at about 
25 mJy, and the strongest high-velocity features 
would only be about 3 mJy (i.e., marginally de¬ 
tectable in a single-epoch GBT spectrum). How¬ 
ever, it is also true that our selection criteria (see 
§3) allowed for the presence of highly asymmetric 
flux ratios in the sample (e.g., an NGC 4258 ana¬ 
logue at a distance of 50 Mpc), yet we found none 
other than NGC 4258 itself. As such, we retain 
the assertion that NGC 4258 is truly anomalous 
in having such a large flux ratio. 


5. Variability 

There are several classes of variability present in 
the megamaser spectra, with different timescales 
and presumed underlying physical causes. We 
qualitatively outline these classes in this section. 

Long-term (^hundreds of days) “bulk variabil¬ 
ity” in the line flux of maser feature sets is seen 
in all sufficiently monitored galaxies. The dynam¬ 
ical timescale for a ~1 pc accretion disk around a 
~10' M 0 black hole is ^10 4 years, so if this bulk 
variability has a dynamical origin, then it likely 
originates from activity much closer to the central 
AGN than any observed masers. Gallimore et al. 
(2001) argue that the megamasers in NGC 1068 
respond to changes in the central power source, 
via a reverberation mechanism. We investigate 
this possibility for several other galaxies in §6. 

Many maser galaxies also display short-term 
(^monthly) flaring variability, where a single 
maser line increases enormously in amplitude, of¬ 
ten by several orders of magnitude over the course 
of only ~a week and lasts for a few weeks. This 
flaring may be caused by the chance alignments of 
individual masing gas clumps in the disk (see, e.g., 
Kartje et al. 1999). In this picture, masing occurs 
in localized clouds which are orbiting ballistically 
in the accretion disk. When one cloud passes in 
front of another while maintaining velocity coher¬ 
ence (as might happen, e.g., for two high-velocity 
clouds on either side of the disk midline), the fore¬ 
ground cloud further amplifies the emission from 
the background cloud, resulting in a rapid increase 
in line luminosity. This provides another potential 
mechanism for the bulk variability, as it could be 
the combined flares of many weak, blended maser 
lines. 

Extremely short-term (intra-day) variability 
that is also uncorrelated among different spec¬ 
tral features has been observed in two megamaser 
galaxies: Circinus (McCallum et al. 2005) and 
NGC 3079 (Vlemmings et al. 2007). This vari¬ 
ability has been attributed to interstellar scintil¬ 
lation, and in §5.2 we present evidence for such 
scintillation in a third megamaser galaxy, ESO 
558-G009. 

We note that our observations are only sen¬ 
sitive to variability on <hourly timescales and 
>monthly timescales. 
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Table 2 



Blue drifts 

Red drifts 

6 p 


Target 

(km s -1 yr -1 ) 

(km s -1 yr -1 ) 

(degrees) 

Reference 


NGC 4258 

-0.140 

± 

0.03 

0.001 

± 

0.004 




Humphreys et al. (2008) 

UGC 3789 

-0.046 

± 

0.04 

0.125 

± 

0.06 




Reid et al. (2013) 

NGC 6264 

0.010 

± 

0.02 

-0.130 

± 

0.01 

< 

0.50 

(b) 

Kuo et al. (2013) 

NGC 6323 

0.030 

± 

0.15 

-0.067 

± 

0.09 

< 

1.50 

(b) 

Kuo et al. (2015) 

Mrk 1419 

0.007 

± 

0.14 

0.052 

± 

0.14 

< 

0.35 

(b) 


NGC 1194 

0.031 

± 

0.13 

0.039 

± 

0.14 

< 

1.55 

(b) 

Litzinger et al. (in prep.) 

NGC 2273 

0.074 

± 

0.23 

-0.011 

± 

0.18 

< 

0.55 

O') 


J0437±2456 

0.036 

± 

0.14 

-0.011 

± 

0.48 

< 

0.55 

0) 


ESO 558-G009 

-0.157 

± 

0.23 

-0.047 

± 

0.22 

< 

6.25 

0) 


IC 2560 

0.011 

± 

0.15 

-0.063 

± 

0.13 

< 

0.55 

(b) 


NGC 5765b 

-0.049 

± 

0.04 

0.008 

± 

0.008 




Gao et al. ( submitted ) 


All -0.036 ±0.014 -0.012 ±0.003 < 0.60 (b) 


Note. —This table lists the mean velocity drifts of high-velocity maser features in the best- 
sampled targets, along with their lcr statistical errors. Values taken from the literature are 
accompanied by the appropriate citations; all other values are new measurements (see §4.2). 
Pitch angles are listed as upper limits, and they are calculated from the velocity drifts of either 
the redshifted (r) or blueshifted (b) features depending on which gives a tighter constraint. 
Values incompatible with the MM98 model have no associated pitch angle. 
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5.1. Dynamic spectra 

One way to effectively visualize both the bulk 
variability and the flaring variability is through 
dynamic spectra. In Figure 5 we present dynamic 
spectra for 9 of our best-sampled sources. To cre¬ 
ate the dynamic spectra, we linearly interpolated 
the flux densities between consecutive GBT spec¬ 
tra, which were taken at a roughly monthly ca¬ 
dence. For the MCP’s monitoring campaign, tar¬ 
gets were not observed during the North American 
summer because atmospheric conditions in Green 
Bank make K-bancl observations inefficient dur¬ 
ing this season. Summer periods with no data are 
blanked. 

Kinematic differences between the systemic and 
high-velocity features, corresponding to differ¬ 
ences in line-of-sight accelerations, are immedi¬ 
ately apparent in the dynamic spectra. Figures 
5(a) and 5(h) match well with Fig. 2 from Braatz 
et al. (2010) and Fig. 1 from Kuo et al. (2013), 
respectively. Further, we note that the systemic 
feature located initially at ~3380 km s -1 in UGC 
3789, which was not used for the distance determi¬ 
nation by Reid et al. (2013) in their acceleration 
analysis for signal-to-noise reasons, shows a clear 
acceleration in the dynamic spectrum. This fea¬ 
ture is offset by about 15 km s -1 from the nearest 
systemic features for which an acceleration was 
measured, so including it would expand the veloc¬ 
ity span of the systemic feature set by ~12% and 
potentially improve the disk model and associated 
distance measurement. 

Along with the kinematic information, the dy¬ 
namic spectra also illustrate how the flux densi¬ 
ties and overall spectral shape change with time. 
If we follow, for instance, the systemic features at 
~3270 km s _1 in UGC 3789, we can see that they 
vary in amplitude by more than an order of mag¬ 
nitude during the ~6-year span of these observa¬ 
tions. We can also see features near this velocity 
appearing and disappearing with time. Several of 
the blueshifted features bracketing 2600 km s _1 , 
on the other hand, remain quite stable in both am¬ 
plitude and structure during the same time range. 
There are also marked differences in feature stabil¬ 
ity among different galaxies; NGC 5765b, for in¬ 
stance, has a very consistent spectrum compared 
to the others. As a result of this spectral stabil¬ 
ity, NGC 5765b has the most precisely-measured 


distance of any MCP galaxy to date (Gao et al. 
submitted). NGC 1194, on the other hand, is ob¬ 
served to be extremely variable; this variability 
has made measurements of this galaxy very chal¬ 
lenging (Litzinger et al. in prep). 

Additionally, we can compare the lifetimes of 
different flaring features in the spectra. The 3270 
km s -1 systemic feature in UGC 3789 flared at 
around day 1700, and it lasted roughly 200 days. 
This duration is considerably longer than that of 
the 3810 km s -1 redslrifted feature, which flared 
around day 1500 but only lasted ~50 days. Com¬ 
pare this to the 1580 km s^ 1 feature in NGC 2273, 
which lasted for at least 400 days, and the 8005 
km s -1 feature in ESO 558-G009, which had a 
duration of ~100 days. 

5.2. Scintillation 

Interstellar scintillation (ISS) in the Galactic 
ionized ISM is considered to be the primary mech¬ 
anism causing the rapid intraday variability ob¬ 
served in pulsars and many extragalactic radio 
sources (predominantly quasars; see, e.g., Bignall 
et al. 2004). For a distant source whose emission is 
undergoing scattering in the turbulent ISM of our 
Galaxy, it is simplest to treat the sum contribution 
from the line-of-sight electron column as originat¬ 
ing from a single thin “scattering screen” located 
a distance D from the Earth. In this picture, tur¬ 
bulence is generated on timescales that are much 
longer than the time it takes a phase-coherent re¬ 
gion of the scattering medium (dubbed a “scintle”) 
to cross the source. That is, the phase variations 
introduced by the screen are essentially “frozen” 
as the screen passes across the line of sight. Thus, 
the scintillation timescale is set by the size and 
transverse velocity of the scattering scintle. 

There are two important ISS regimes separated 
by a “transition frequency” v t \ the weak (v > v t ) 
and strong (v < v t ) scattering limits. We give a 
brief overview of some relevant properties of these 
limits here; for a review of this topic, see Narayan 
(1992) and references therein. 

In the weak scattering limit, the size of the scin¬ 
tle is of order the Fresnel scale, defined to be the 
transverse distance from the line of sight to a point 
through which the increase in path length from 
the source to observer (compared to the direct, 
line-of-sight path) results in a phase change of 1 
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Fig. 5.— Dynamic spectra for our best-sampled disk megamasers. For ease of viewing, the three sets of 
features have been split up and the spectral regions in between (which are devoid of maser features) are not 
shown. The color scale maps to the logarithm of the flux density, as shown in the colorbar on the right. 
Individual observation dates are indicated by white tick marks near the bottom of each plot, and day zero 
is set as the date of the first observation (see Figure 1). Velocities are measured in the heliocentric frame, 
using the optical velocity convention. 
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Fig. 5.— (continued) 
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Fig. 5.— (continued) 
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Fig. 5.— (continued) 
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radian. For a source at infinity and an observ¬ 
ing wavelength A <C D, the Fresnel scale is given 
by r F = yOcD/2 7 r. If the scattering screen has 
transverse velocity (relative to the Earth) of v, the 
variability timescale will be r « rp/v. 

In the strong scattering limit, the scintle has 
a characteristic size called the diffractive scale, 
7"diff. This length scale functions equivalently to 
the Fresnel scale in weak scintillation (i.e., the 
RMS phase difference between two points on the 
screen separated by a distance rdiff is approxi¬ 
mately 1 radian), but the physical origin of the size 
scale is different. In the strong scattering regime, 
the value of rdiff is determined by the turbulent 
properties of the ISM plasma rather than by the 
geometry of the observer-screen-source setup. We 
thus have rdiff r F for strong scattering, while 
r F -C rdiff- for weak scattering. The scintillation 
timescale will then be r k i'd\s/v. The strong 
scattering regime can be further subdivided into 
two different types of strong scattering, diffractive 
and refractive. Refractive scintillation occurs on 
much longer timescales (^days) than diffractive 
scintillation, so it is not relevant for this study. 

A standard measure of variability strength is 
the modulation index, /i = a/(S), where er is 
the standard deviation of the observed amplitude 
and ( S ) is its average value. The modulation 
index for a point source undergoing weak scat¬ 
tering is roughly the ratio of the Fresnel to the 
diffractive scale, n « (r F /rdiff) 5 ^ 6 (Narayan 1992). 
For diffractive scintillation, the modulation index 
should be unity. In the case of an extended source 
(i.e., a source with an angular size larger than 
the diffractive scale), the diffractive scintillation 
is said to be “quenched,” since the resolved source 
is effectively diluting the variability amplitude by 
averaging the phase fluctuations over several adja¬ 
cent scintles. An extended source of angular size 
8 will have a modulation index given by d^g/d. 

ISS has been proposed as an explanation for 
the extremely rapid (intra-hour) variability ob¬ 
served in the 22 GHz maser spectra from the Circi- 
nus galaxy and NGC 3079. Vlennnings et al. 
(2007) use the high Galactic latitude of NGC 
3079 (b = +48.36°) to justify their assumption 
of weak scintillation. From a measured character¬ 
istic timescale of r « 1000 s, corresponding to the 
crossing time for the Fresnel scale, they calculate 
a distance to the scattering screen of D « 25 pc. 


McCallum et al. (2005) measured the timescale 
in Circinus to be t « 700 s, but were unable to 
say definitively whether the variability was caused 
by weak scintillation in a nearby screen (D « 20 
pc) or quenched diffractive scintillation in a more 
distant screen (D ss 230-1000 pc). Followup ob¬ 
servations from McCallum et al. (2007) showed 
spectral variations that lent strong support to the 
diffractive scintillation interpretation, and they 
further uncovered longer-timescale (~1 day) vari¬ 
ations consistent with refractive scintillation. 

5.2.1. Scintillation in ESO 558-G009 

We present here observations of the third mega¬ 
maser galaxy observed to show signs of ISS. Figure 
6 shows light curves for two epochs of the galaxy 
ESO 558-G009 on which we’ve applied our scintil¬ 
lation analysis. These epochs were chosen because 
of their long observation durations (> 3 hours 
each) and because they both contained the same 
strong systemic maser feature (> 0.15 Jy), which 
was detectable in a single 5-minute scan. We ex¬ 
amined the spectra for all the other megamasers 
that met these same criteria (long-duration obser¬ 
vation and strong maser feature), but only ESO 
558-G009 showed significant variability. Figure 6 
also shows the discrete autocorrelation functions 
(DACFs) for both of the light curves, calculated 
using the technique outlined by Edelson & Krolik 
(1988). The dates of the observations are listed in 
Table 3. 

The light curves show variability timescales on 
the order of ~2100 s, during which the peak flux 
can vary by a factor of ~3; this is comparable to 
the amplitude modulations observed in the quasar 
J1819+3845, the extragalactic source exhibiting 
the strongest ISS-induced continuum variability 
(Dennett-Thorpe & de Bruyn 2002). Though it’s 
possible for pointing errors or a changing atmo¬ 
spheric opacity to cause the amplitudes of spectral 
features to vary with time, we don’t expect these 
effects to exceed ~20%. Further, if such factors 
were the cause of the observed amplitude changes 
then we would expect to see them across spectra 
of all galaxies, which is not the case. As a final 
check, we measured the total flux of the systemic 
and high-velocity features (outside of the targeted 
line) over time during the observations, and we 
found that it is constant to within ~15% through¬ 
out a single observing session. 
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If the variability were intrinsic to the maser, 
then such large amplitude changes must result 
from increases in the maser gain path that are of 
order the gain length, £, which for an unsaturated 
maser is the path length corresponding to an e- 
fold increase in amplification (i.e., it is the length 
over which the optical depth r changes by ~1). 
For conditions typical of those found in megamaser 
disks, £ 1 AU (Greenhill et al. 1997a). Given 

the light-travel distance of ~4 AU derived from the 
characteristic timescale, the observed variability 
would require changes in the gain path to propa¬ 
gate at approximately the speed of light. Barring 
radiative pumping (which is not expected to be 
important in these systems; see, e.g., Lo 2005), we 
do not know of any mechanism capable of driving 
such rapid changes. This leaves foreground scin¬ 
tillation as the best available explanation for the 
variability. 

Following Rickett et al. (2002), we define the 
characteristic observed scintillation timescale, r, 
to be the half-width at half-maximum (HWHM) 
of the autocorrelation function. If the masers be¬ 
have as a point source (i.e., if their angular size 
is smaller than the angular size of the scattering 
screen, 6 < 6 S ), then a measurement of r allows us 
to establish a characteristic size, r s , for the scat¬ 
tering screen (i.e., the size of a scintle) of 

r s = v s t. (10) 

Here, v s is the transverse velocity of the screen rel¬ 
ative to Earth. In Appendix B we have outlined 
how this transverse velocity is obtained for an in¬ 
dividual observation, using a model that combines 
the Earth’s orbital motion and the Sun’s peculiar 
and orbital motion. Table 3 lists v s for each ob¬ 
servation, assuming a nearby (D < 100 pc) screen; 
the measured values for r are also listed. 

Our model assumes that the scattering screen 
itself has no peculiar motion. From time-delay 
measurements of the intra-day variability in the 
quasar J1819+3845, Dennett-Thorpe & de Bruyn 
(2002) found that the scattering screen (for that 
target) must have a transverse peculiar velocity of 
about 25 km s _1 . McCallum et al. (2009) used the 
same technique to place a lower limit of 22 km s -1 
on the transverse velocity of the ISM along the line 
of sight to Circinus. We have no reason to expect 
that the scattering screen towards ESO 558-G009 


should behave any differently. However, with two 
free parameters already in the model (D and r s ) 
and only two measurements, we have no room to 
add the two additional parameters that would be 
necessary to properly account for peculiar motion. 
We are thus only able to place relatively crude 
constraints on the model parameters. Figure 7 
shows these constraints, with the more relevant 
9 S — Ts/D plotted in place of r s . We can see that 
our measurements, which have a formal “best fit” 
at about D ss 70 pc and 9 S ss 5 ^.as, are compatible 
with a wide range of parameters. A scintle angular 
size of 5 /ias corresponds to a lower limit on the 
maser brightness temperature of ~3 x 10 13 K. 

If we assume that the scintillation occurs in the 
weak scattering regime, then we have r s ss tf, and 
we can use the Fresnel scale to determine D. Do¬ 
ing so yields a distance to the scattering screen be¬ 
tween 40 and 50 pc. From Walker (1998), we can 
use the modulation index to determine the tran¬ 
sition frequency. Between the two observations, 
/z « 0.5, so we obtain v t ~ 13.6 GHz. 

Like Circinus, ESO 558-G009 is located near 
the plane of the Galaxy (b = —6.96°), so we 
would expect to see greater-than-average scatter¬ 
ing along this line of sight. From the NE2001 
model for the electron density along different lines 
of sight in the Milky Way (Cordes & Lazio 2002), 
the transition frequency between weak and strong 
scintillation towards ESO 558-G009 should actu¬ 
ally be about 30 GHz; since this is higher than the 
observing frequency of 22 GHz, it would put us in 
the strong limit. We note that the Cordes & Lazio 
(2002) model attempts to map the Galactic elec¬ 
tron density in a primarily spatially smooth man¬ 
ner, while the true distribution is known to have 
mesoscale and microscale structure. As such, we 
expect significant model uncertainties along any 
specific line of sight. 

In the strong scintillation regime the measured 
timescale maps to the angular size of the source 
rather than to that of a scintle. The modulation 
index should be equal to the ratio 9dis/9 Sl so a 
modulation index of /i ss 0.5 (see Table 3) indi¬ 
cates that the angular size of the maser must be 
a factor of ~2 larger than that of the scintle. For 
a screen distance of 70 pc we have 6 S « 5 /ias. 
For the ESO 558-G009 distance of 110 Mpc, we 
thus obtain an approximate physical size of the 
masing region of ~1100 AU. This is comparable 
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Fig. 6. — Light curves (left) and discrete autocorrelation functions (DACF, right) for two observations of 
ESO 558-G009. In the light curves, the LCP and RCP peak flux densities of the 7590 km s -1 maser line are 
plotted (with circles and squares, respectively) at the ~ 5-minute cadence corresponding to individual nod 
scan pairs. The dotted vertical line in the DACF marks the location of r (i.e., where the DACF drops to a 
value of 0.5). 
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to the 0.001-0.006 pc clump sizes estimated by 
Kondratko et al. (2005) for the disk of NGC 3079. 

6. Testing for disk reverberation 

Claussen & Lo (1986) noted that the apparent 
systematic flux variations in the nuclear masers 
in NGC 1068 suggested that the masers share a 
common pumping source. If all masers in a given 
galaxy are powered by a common source, presum¬ 
ably at the nucleus, then we would expect vari¬ 
ability in the power source to propagate to the 
maser system. This variation would reverberate 
through the masing disk at some propagation ve¬ 
locity which, if it is on the order of the speed of 
light, would be fast enough to pass through the 
entire masing portion of the disk on timescales 
of a year or two. Gallimore et al. (2001) mea¬ 
sured a correlation between the variability of red- 
shifted and blueshifted maser features in NGC 
1068, which they used to argue that the masers 
respond to variability in the central engine. 

Since the fiducial picture of circumnuclear 
megamaser disk geometry (for a Keplerian ro¬ 
tation curve) allows us to uniquely associate any 
high-velocity maser feature with a radial location 
within the accretion disk, we attempt here to de¬ 
tect the propagation of some signal through the 
masing disks of our best-sampled targets. A mea¬ 
surement of disk reverberation would not only lend 
support to the idea of a common pumping source, 
but it could also potentially enable an independent 
means of measuring the mass of the central SMBH 
and the distance to the host galaxy (provided the 
propagation velocity of the reverberation signal is 
known). If we denote the outward propagation 
speed as v s , then a reverberation signal passing 
through the spectrum at a rate v corresponds to 
a black hole mass of 




v s (v - v 0 ) 3 
2 Gv 


( 11 ) 


Here, v is the observed velocity (i.e., as seen in the 
spectrum) and vq is the velocity of the dynamic 
center (i.e., the motion of the black hole itself, 
which is presumably almost identical to the reces¬ 
sion velocity of the galaxy). We note that v will in 
general be a function of v; that is, for a constant 
value of v s the rate at which the reverberation sig¬ 
nal passes through the spectrum depends on where 


in the spectrum it is located. Once the black hole 
mass is known, the distance to the galaxy can be 
determined by comparing the angular orbital radii 
of the maser spots (measured using VLBI) to the 
orbital radii calculated using the single-dish spec¬ 
tra (from r = GMbb/v 2 ). 

6.1. Extracting a reverberation signal 

Here we outline the procedure used to check for 
the spectral signature of radially-propagating exci¬ 
tation in a time series of GBT disk maser spectra. 
The relevant parameters are the mass of the cen¬ 
tral black hole, Mbh, the recession velocity of the 
dynamic center, Vq, and the propagation speed of 
the signal, v s . The observed response of a high- 
velocity maser offset by a distance D (see bottom 
panel in Figure 8) is delayed by D/v s relative to 
the response of all systemic masers. 

We subtract a weighted average spectrum (see 
§3) of the target from each epoch to remove stable 
(i.e., non-propagating) high-velocity features from 
each spectrum. We then map each velocity chan¬ 
nel, Vi, to a radial position, r.j, within the maser 
disk. The mapping assumes that the high-velocity 
maser spots are all located on the midline of the 
disk, and that they are all on circular Keplerian 
orbits: 


GMbb 
(Vi - v 0 ) 2 ' 


( 12 ) 


We refer to the original GBT spectra as the “veloc¬ 
ity spectra” and the new, radially-mapped spectra 
as the “radial spectra.” An example of these two 
for the source UGC 3789 is shown in Figure 8. 

To account for the time delay between the 
detection of a propagating signal in consecutive 
epochs, each radial spectrum is temporally shifted 
according to the signal propagation speed and that 
spectrum’s date of observation, relative to some 
reference epoch. For simplicity, we have defined 
the temporal zeropoint to be the date of the first 
observation, given in Figure 1. This process is il¬ 
lustrated in Figure 9. 

After shifting, the radial spectra are then aver¬ 
aged over all epochs. If a target has been observed 
for N epochs, each of which has an associated ra¬ 
dial spectrum S n {r,t n ), then this procedure can 
be written as 
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Table 3 


Date 

V 

(km s -1 ) 

v s 

(km s -1 ) 

(5) 

(Jy) 

a 

(Jy) 


r 

(hours) 

2011 Oct 09 

7589.8 

24.4 

0.186 

0.082 

0.44 

0.60 ±0.07 

2012 Feb 21 

7590.3 

27.4 

0.166 

0.088 

0.53 

0.58 ±0.08 


Note. —Scintillation parameters for ESO 558-G009. The column titled v 
lists the Doppler velocity for the targeted line, v s is the modeled transverse 
velocity at the observation date (see Appendix B for details), (S) is the 
mean flux density for the line during the observation, a is its standard 
deviation, fj, is the modulation index, and r is the measured characteristic 
variability time. 



Fig. 7. — Constraints on the angular size of the scintles and the distance to the scattering screen along the line 
of sight to ESO 558-G009. The solid line (with 3<r error in blue) shows the constraint from the 2011 October 
9 observation, and the dashed line (with 3cr error in red) shows the constraint from the 2012 February 21 
observation. The plotted errors account only for the statistical errors arising from the measurement of r; no 
systematic errors from the velocity modeling are included. 
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Fig. 8.- Illustration of the conversion between a velocity spectrum (top) and a radial spectrum (bottom), 
using Equation 12. The dashed line in the upper spectrum shows the recession velocity of the system, and the 
black point in the lower spectrum shows the location of the SMBH. The blueshifted portion of each spectrum 
is plotted in blue, while the redshifted portion is plotted in red. The source chosen for this example is the 
galaxy UGC 3789. 
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Fig. 9.— These plots show the radial spectra before (left) and after (right) accounting for the time delay 
caused by the propagation of the signal. The black spectra are real spectra of UGC 3789, and the red line 
includes the artificially injected 10 mJy signal. In the panel on the left, we can see that the artificial signal 
is propagating outwards with time. In the panel on the right, the spectra have been temporally shifted using 
v s = c; as a result, when stacking these spectra the signal will add coherently. In both panels, the spectra 
have been vertically offset by an amount proportional to the time between observations; the time since the 
first observation is shown on the right axis. The radial zeropoint corresponds to the position of the SMBH 
at the time of the first observation. Only the redshifted high-velocity features are shown in these plots. 
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S(r) = jy S n (r - v s t n ,tn)- (13) 

n—1 

Here, S(r) is the final combined radial spectrum. 
The radial zeropoint is defined to be the center of 
the disk (i.e., the location of the SMBH) at the 
date of the first observation. 

The purpose of this procedure is to stack spec¬ 
tra in such a way that a radially-propagating sig¬ 
nal will add coherently across all epochs. 

6.2. Sensitivity 

The sensitivity of our method depends on sev¬ 
eral factors, including the number of epochs and 
overall time baseline of observation, as well as the 
intrinsic variability of the target. We restrict our 
analysis to well-sampled (i.e., >20 epochs of obser¬ 
vation) galaxies that have reliably measured black 
hole masses (see Kuo et al. 2011). 

Some of these targets are more variable than 
others. In general, the more flaring a source dis¬ 
plays, the less sensitive this measurement will be. 
Flaring events are not removed well when sub¬ 
tracting an epoch-averaged spectrum, and so suf¬ 
ficiently strong flares can appear as false positives 
in the final radial spectrum. 

Although we’ve chosen to test only those 
sources for which Mbh is known to ~10% or bet¬ 
ter, it’s possible that our method requires the 
value to be even more precisely known to ensure 
recovery of a propagating signal. To test the sensi¬ 
tivity of our method on the input values of v s and 
Mbh, we injected an artificial propagating signal 
into a series of spectra. The signal was a Gaussian 
pulse of fixed width and amplitude, propagating 
with a fixed velocity from a black hole of known 
mass. We found that the tolerance threshold for 
both v s and Mbh was approximately 5%; if either 
of these inputs is off from the true value by more 
than this amount, the signal is not recovered or is 
severely degraded. 

6.3. Discussion 

We tested for reverberation in the six maser 
galaxies listed in Table 4. For each galaxy, we 
checked for signals propagating at velocity in¬ 
crements of 0.01c, with minimum and maximum 


propagation velocities of 0.8c and 1.2c, respec¬ 
tively 1 . We also adjusted the black hole masses 
within a range ±20% of the measured value, in 
increments of 1%. No reverberation signals were 
detected in any of the galaxies, with limiting flux 
densities listed in Table 4. Given that the spec¬ 
tra for these galaxies typically vary at the Mens 
of mJy level (see §5), we can see that any contri¬ 
bution from a propagating signal must constitute 
only a small (< 10%) fraction of the total variabil¬ 
ity. 

The detection thresholds listed in Table 4 are 
simply the 3er noise levels in the final combined 
spectra. We emphasize that this threshold gives 
only a limiting value for a signal that is perpet¬ 
ually coherent (i.e., always maintains its profile 
shape and moves at constant velocity) and that 
is present in all available spectra (i.e., it does not 
fade in and out as it propagates). This procedure 
is less sensitive to a more complex signal. 

Furthermore, we note that the timescale for 
variability of the pumping source influences our 
measurements. If the source doesn’t vary much 
over the Mew-year timescales probed by these 
data, then the signal won’t be radially localized 
and our technique will not help to detect it. 

7. Magnetic field strengths from Zeeman 
splitting 

Magnetic fields in AGN accretion disks are 
thought to drive several important physical pro¬ 
cesses. The magnetorotational instability (MRI), 
first described in a general astrophysical context 
by Balbus & Hawley (1991), is likely the primary 
means by which angular momentum is transported 
in accretion disks. Magnetic fields are also neces¬ 
sary for launching outflows, from the classic MHD 
disk wind (Blandford & Payne 1982) to more mod- 


1 We actually investigated propagation speeds down to 0.0c, 
but the sensitivity of the method starts to drop consider¬ 
ably below a certain speed. This is because the individual 
maser features — which in general aren’t perfectly matched 
to the average spectrum, so they don’t subtract out well 
- begin to add semi-coherently, rather than averaging out 
like noise. To give an example, the 3cr threshold for UGC 
3789, which is about 0.8 mJy for propagation speeds be¬ 
tween 0.8c and 1.2c, increases to ~5 mJy for a propagation 
speed of 0.5c. This also makes it more difficult to differen¬ 
tiate between a propagating signal and a coherently-added 
maser feature, so we only quote sensitivities between 0.8c 
and 1.2c for Table 4. 
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Table 4 


Target 

Vq (km s : ) 

M bh (10 7 M 0 ) 

Epochs 

Threshold (mJy) 

UGC 3789 

3262 

1.04 

58 

0.8 

Mrk 1419 

4954 

1.16 

55 

1.4 

NGC 6323 

7829 

0.94 

44 

1.2 

NGC 1194 

4063 

6.5 

43 

4.1 

NGC 2273 

1832 

0.75 

38 

1.8 

NGC 6264 

10194 

2.91 

28 

0.8 


Note. —Galaxies tested for a reverberation signal. The threshold column 
lists the 3cr detection cutoffs; a signal stronger than this value would be 
classified as a detection. Note that the velocity of the dynamic center (v 0 ) 
need not be the same as the recession velocity of the galaxy listed in Table 
1, as the Vo values were obtained by fitting Keplerian rotation curves to 
position-velocity data (Kuo et al. 2011). 


ern incarnations that also incorporate radiation 
pressure (e.g., Keating et al. 2012). In this sec¬ 
tion, we use measurements of the Zeeman effect 
to place limits on the magnetic field strength in 
several megamaser disks. 

The maser emission that we observe at 22 GHz 
arises from one or more of the six hyperfine transi¬ 
tions of the 6i6-523 rotational transition of the wa¬ 
ter molecule (see Fiebig & Guesten 1989). Since 
this molecule is non-paramagnetic, Zeeman split¬ 
ting of these hyperfine energy levels arises from 
the coupling between the nuclear magnetic mo¬ 
ments and an external magnetic held. This causes 
the effect to be much weaker (by a factor of ~10 3 ) 
in water than in molecules such as OH, where the 
unpaired electron’s spin couples with the magnetic 
held. The drastic difference in magnitude arises 
because the Bohr magneton and the nuclear mag¬ 
neton differ by the ratio of the electron to the nu¬ 
cleon mass, m e /m p « 1/1836 (Heiles et al. 1993). 

An external magnetic held causes each hyper- 
hne level to split into three groups of lines: the 
7 r components and the cr^ components, corre¬ 
sponding to magnetic quantum number changes 
of AMp = 0 and AM j? = ±1, respectively (Mod- 
jaz et al. 2005). The a ± components are circu¬ 
larly polarized about the magnetic held direction, 
and they are symmetrically offset from the parent 
frequency. For weak magnetic helds (i.e., B < 1 


Gauss), this frequency offset is small compared to 
the line width; typically (Av z / Avl) ~ 10 -3 -10 -4 . 

7.1. Method 

In principle, the measured frequency difference 
between the left and right circular polarizations 
(corresponding to cr + and a ~, respectively) allows 
us to determine the line-of-sight component of the 
magnetic held at the location of the maser spot. 
Since the offset is small compared to the width of 
the line profile, the Stokes V profile (given by V = 
[LCP—RCP]/2) is proportional to the derivative of 
the Stokes / profile (given by I = [LCP+RCP]/2). 
This leads to a characteristic S-shape of the Stokes 
V prohle (see, e.g., Vlemmings et al. 2001, Fig. 2). 

Modjaz et al. (2005) conducted a series of 
Monte Carlo simulations which established that 
the RMS sensitivity to the line-of-sight compo¬ 
nent of the magnetic held from a single maser line 
is consistent with what one would expect from 
a statistical treatment (see, e.g., Lenz & Ayres 
1992), namely: 


<?b 


Av l 
2 A 



(14) 


Here, Avl is the FWHM line width, S/N is the 
Stokes I signal-to-noise ratio, and A is the Zeeman 
splitting coefficient (which is different for each hy- 
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perfine transition). After numerically solving the 
radiative transfer and rate equations for magne¬ 
tized water masers, Nedoluha & Watson (1992) 
found that a value for A of 0.020 km s _1 G -1 
was most appropriate for the merging of the three 
dominant hyperfine components. This value as¬ 
sumes that the three strongest hyperfine lines all 
contribute to a given observed maser line, and de¬ 
viations from this value never exceeded a factor 
of ~2 across the range of parameter space inves¬ 
tigated in Nedoluha & Watson (1992). We thus 
adopt A = 0.02 km s _1 G -1 for our calculations, 
which in general follow the same procedure out¬ 
lined in Modjaz et al. (2005). 

For extragalactic sources, only three efforts to 
measure magnetic held strengths using the Zee- 
man effect in H 2 O megamasers have been pub¬ 
lished. Modjaz et al. (2005) placed a lcr up¬ 
per limit of 30 mG on the radial component of 
the magnetic held in NGC 4258, using a cross¬ 
correlation method to handle the heavy blending 
of the spectral features. Vlenmiings et al. (2007) 
used the same technique on NGC 3079, obtain¬ 
ing an upper limit of 11 mG for the blueshifted 
features. Both studies also measured limits for 
strong, isolated maser components, and combined 
these results with those from the cross-correlation 
method. Additionally, McCallum et al. (2007) 
measured isolated lines to place a lcr upper limit of 
50 mG on the toroidal component of the magnetic 
held in the Circinus galaxy. 

7.2. Measurements 

The most sensitive test for Zeeman splitting us¬ 
ing individual (i.e., non-blended) maser lines oc¬ 
curs on lines that are both strong (large signal-to- 
noise) and narrow (small Avl)- We therefore fo¬ 
cused our test on strong ( S/N > 50) haring events. 

For each selected maser hare, we separately re¬ 
duced the LCP and RCP spectra without applying 
Hanning smoothing; this process retains the full 
spectral resolution. To compensate for errors in 
flux scale calibration, the peak value of the RCP 
spectrum was scaled to the value of the LCP spec¬ 
trum prior to computing either the Stokes I or 
Stokes V spectra. Typical scaling offsets were of 
order 10%. We note that the absolute intensity 
scale is unimportant for these measurements. 

We did not detect Zeeman splitting in any of 


the maser lines, so our results here yield only up¬ 
per limits on the magnetic held strengths. These 
results are summarized in Table 5, and an example 
measurement (from NGC 1194) is shown in Figure 
10 . 

7.3. Discussion 

Since the Zeeman measurements are only sen¬ 
sitive to the line-of-sight magnetic held, B«, the 
high-velocity and systemic lines measure different 
equatorial components of this held. The high- 
velocity features measure the toroidal component 
of the held, B toI , while the systemic features mea¬ 
sure the radial component, B r . None of the fea¬ 
tures directly measure the poloidal component of 
the magnetic held, but an appropriate model (see, 
e.g., Hawley et al. 1996) can estimate its magni¬ 
tude using the values of the toroidal and radial 
components. 

Even without knowledge of the poloidal com¬ 
ponent, we can still use the derived upper limits 
to constrain the support mechanism for the accre¬ 
tion disks. This is because only the components 
of the held that thread through the disk (i.e., only 
the radial and toroidal components) can provide 
vertical pressure support. For typical maser con¬ 
ditions of n ~ 10 9 cm -3 and T ~ 1000 K, the gas 
pressure amounts to roughly 10 -4 erg cm -3 . The 
equivalent support from magnetic pressure would 
require a ~50 mG magnetic held, which is com¬ 
parable to (though still slightly below) our most 
stringent limits. It is worth noting that these num¬ 
bers are also comparable to the ~100 mG upper 
limit imposed by hydrostatic equilibrium for the 
disk thickness measured by Argon et al. (2007) in 
NGC 4258. 

8. Summary 

We have addressed several new scientific ques¬ 
tions that can be explored using the MCP’s exten¬ 
sive monitoring campaign of 22 GHz disk mega¬ 
maser spectra with the GET. The spectra in this 
dataset are unique in their ability to probe the ac¬ 
cretion disks of nearby AGN at sub-parsec scales, 
and the dataset itself is unmatched in the sensi¬ 
tivity and time coverage for each target. 

1. We present a comprehensive collection of 
Keplerian disk megamaser spectra. We also 
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Velocity (km/s) 


Fig. 10. - GBT spectrum of NGC 1194, taken on 2011 December 30. Inset are the Stokes I and V profiles 
for the 4097.2 km s -1 feature (left) and the 4751.7 km s -1 feature (right). The black dashed lines in the 
Stokes V plots show the lcr RMS level for this spectrum. No Zeeman profile is evident for either of these 
lines; limits are given in Table 5. 











Table 5 


Target 

Date 

-Mbh 
(10 7 M 0 ) 

V 

(km s -1 ) 

Vrot 

(km s -1 ) 

Peak 

(mjy) 

S/N 

Av l 

(km s -1 ) 

B \\ 

(inG) 

Radius 

(pc) 

NGC 1194 

2007 Dec 26 

6.5“ 

4757.6 

694.6 

340 

141 

0.58 

<100 (t) 

0.58 

NGC 1194 

2010 Apr 10 

6.5 

4146.4 

83.4 

210 

97 

0.87 

<220 (r) 

- 

NGC 1194 

2011 Dec 30 

6.5 

4097.2 

34.2 

1020 

330 

0.96 

<73 (r) 

- 

NGC 1194 

2011 Dec 30 

6.5 

4751.7 

688.7 

800 

259 

0.95 

<91 (t) 

0.59 

NGC 2273 

2009 Dec 12 

0.75“ 

1582.5 

-249.5 

240 

115 

0.72 

<160 (t) 

0.52 

NGC 3393 

2006 Apr 28 

3.1 6 

4050.9 

300.9 

230 

95 

0.84 

<220 (t) 

1.48 

NGC 3393 

2006 Dec 6 

3.1 

4260.8 

510.8 

350 

94 

1.1 

<300 (t) 

0.51 

UGC 3789 

2010 Dec 20 

1.04“ 

3273.0 

11.0 

190 

75 

0.74 

<250 (r) 

- 

NGC 6323 

2008 Mar 25 

0.94“ 

7395.2 

-433.8 

180 

86 

1.0 

<300 (t) 

0.21 

ESO 558-G009 

2013 Apr 22 

1.8“ 

8003.7 

329.7 

490 

81 

0.99 

<310 (t) 

0.71 

Mrk 1419 

2007 Apr 14 

1.16“ 

5330.8 

376.8 

220 

56 

1.6 

<720 (t) 

0.35 


Note.— Maser lines tested for Zeeman splitting. For flaring lines appearing in more than one epoch, the listed observation 
date is that which yields the best upper limit on the line-of-sight component of the magnetic field. In addition to the Doppler 
velocity (a), we list the rotation velocity (foot = v — vo', blueshifted lines are negative, vo is the velocity of the dynamic 
center) and the measured line width (A»i) for each line. For all lines, By is quoted as a la upper limit, and the letters in 
parentheses indicate whether the measurement is sensitive to the toroidal (t) component or the radial (r) component. For 
limits on toroidal magnetic field components, the radius column gives the corresponding radial location in the disk at which 
the limit holds. 

“Kuo et al. (2011) 

^Kondratko et al. (2008) 

“Gao et al. ( submitted ) 


present dynamic spectra for the most heavily 
monitored of these sources. 

2. We find that the redshifted high-velocity 
maser features are brighter, on average, than 
the blueshifted features for our sample of 32 
megamaser disks. This asymmetry is pre¬ 
dicted by the spiral shock model of MM98. 

3. We also test the MM98 prediction that the 
high-velocity features should exhibit nonzero 
line-of-sight velocity drifts. We find no sys¬ 
tematic drifts. Furthermore, the statistically 
significant detection of both positive and 
negative velocity drifts within the same set 
of features (as we have for several sources) 
is inconsistent with the MM98 model’s pre¬ 
dictions. 

4. We argue that the intra-day variability 
observed in ESO 558-G009 is most likely 
caused by ISS, and we derive parameters of 
the scattering screen under different assump¬ 
tions about the scattering regime. Though 
the measurements are currently sparse, we 


find that they are most consistent with a rel¬ 
atively nearby (~70 pc) scattering screen. 

5. We test six maser systems for a radially- 
propagating change in maser activity, which 
could be the result of variable output from 
the central engine. No such signal is detected 
in any of the galaxies. 

6. We measure upper limits on the toroidal and 
radial magnetic field strengths in the accre¬ 
tion disks of 7 galaxies using the Zeeman ef¬ 
fect, and we find that the magnetic fields 
must be less than several hundred rnG in 
each case. This is beginning to probe the 
regime where the magnetic pressure becomes 
comparable to the gas pressure in the disk. 

We acknowledge Fred Schwab for guidance on 
the statistical analysis, Ken Kellerman for a dis¬ 
cussion of variability timescales, and Scott Suri- 
ano for helpful clarifications regarding disk wind 
phenomenology. The National Radio Astronomy 
Observatory is a facility of the National Science 
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Foundation operated under cooperative agree¬ 
ment by Associated Universities, Inc. This re¬ 
search was supported in part by an ARCS/MWC 
Scholar Award. This research has made use of 
the NASA/IPAC Extragalactic Database (NED), 
which is operated by the Jet Propulsion Labo¬ 
ratory, California Institute of Technology, under 
contract with the National Aeronautics and Space 
Administration. 
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A. Kernel density estimation (KDE) 

In essence, the KDE technique as used in this paper is simply an alternative to a traditional histogram 
(though in each case we have shown it alongside such a histogram). In a standard histogram, a single data 
point falls into a “bin” of width h, unit height, and fixed edgepoints. The bin width is usually determined 
by the sample size and spread, with the optimal result being a compromise between data resolution and 
population per bin. The bin edgepoints, however, are often more arbitrarily defined. The KDE approach 
solves this issue by eliminating the use of bins; instead, each data point is represented by a “kernel” of some 
predefined functional form. In Figures 2 and 3 we used a Gaussian kernel of the form 

K(u) = ~^e~ u2/2 . (Al) 

This kernel has been scaled by h relative to a normal (i.e., unit area) Gaussian kernel, such that the area 
under the curve for a given data point matches what would be found in a typical histogram of bin width h. 
The final kernel density estimator is then just a sum of the kernels for all data points, which can be written 
as 


= f>'■ < A2 > 

Here, N is the number of data points, A, is the center of the kernel for data point i (i.e., the value of that 
data point), and h is the kernel width. In our case, X, is the isotropic luminosity (for Figure 2) or the 
logarithm of the flux ratio (for Figure 3) for a single source. Under the assumption that our underlying 
distribution is at least approximately Gaussian, we have used the bin width derived by Silverman (page 45, 
equation 3.28): 



Here, a is the standard deviation of the sample. 


(A3) 


B. Transverse motion along the line of sight to ESO 558-G009 

Our goal is to transform from the Galactic Cartesian coordinate system (A, Y, Z) to a coordinate system 
(x, y, z) where the line of sight to ESO 558-G009 is aligned with the 2 -axis. We define 2 : to be pointing away 
from ESO 558-G009 and y to be the projection of the North Ecliptic Pole onto the plane perpendicular to 
z. The unit vector x is then defined to be x = z x y. 

The standard spherical Galactic coordinates (£, b) can be converted to Galactic Cartesian unit vectors 
using the transformation: 


X = cos (£) cos (b) (Bl) 

Y = sin(f) cos (b) 

Z = sin(6) 

We can thus define a unit vector r = (A, Y. Z) that points in the direction of any Galactic coordinate location 
(£,b). 

The North Ecliptic Pole has Galactic coordinates (£, b) = (96.3840,29.8117), with corresponding unit 
vector f NEP = (-0.0965,0.8623,0.4972). The coordinates for ESO 558-G009 are (£, b) = (233.6609, -6.9598), 
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with unit vector teso = (—0.5882, —0.7996, —0.1212). From our description above of the desired coordinate 
system, we have the following expressions for the coordinate unit vectors: 


r N EP x r E so 
FNEP X r E so| 
y = T’ESO X X 
z = -f'ESO 

These evaluate to x = (0.4064, -0.4218,0.8105), y = (-0.6992,0.4275,0.5731), and z = (0.5882,0.7996,0.1212). 
We’ll henceforth refer to this new coordinate system as the “source” coordinate system. 

B.l. Solar motion with respect to the LSR 

The first component of the transverse motion comes from the Sun’s deviation from its orbital mo¬ 
tion. From Co§kunoglu et al. (2011), the Sun’s peculiar motion relative to the LSR has components 
v © = (8.50,13.38,6.49) km s” 1 , with magnitude v© = 17.13 km s -1 and corresponding unit vector 
f© = (0.4962,0.7811,0.3789). The parallel and perpendicular components of this velocity are then sim¬ 
ply its projections onto the coordinate axes: 


%,|| = v e ' z (B3) 

v Q ,± = \J(v e • xf + (u© • ijf 

These evaluate to u© m = 16.48 km s _1 and u©,© = 4.66 km s _1 , with source components (v x ,v y ,v z )@ = 
(3.07,3.50,16.48) kin s" 1 . 

B.2. Earth’s orbital motion 

The second component of the transverse motion comes from the Earth’s orbit around the Sun. For 
simplicity, we'll model this orbit as circular about the North Ecliptic Pole, with orbital velocity v© = 30 
km s _1 . If we define a position angle <j> = u>t measured clockwise from the negative rr-axis, then we can 
decompose the Earth’s orbital motion into the source components: 


v x , ®(t) = v e sin((/> 0 + ojt) (B4) 

v y, ®{t) = f© cos(0 o + wt) cos (i) 

v z, ®(t) = - v © cos(<^o + wt) sin(i) 

Here, u> is the orbital angular frequency of the Earth, i = 7r/2 — cos -1 (y ■ tnep) is the inclination of the 
orbit relative to the line of sight to ESO 558-G009 (in our case, i = 46.1°), and <j >o is an initial position angle 
that must be calibrated based on the known motion of the Earth. 

On the vernal equinox (the origin of the ecliptic longitude), the Earth is moving towards ecliptic co¬ 
ordinates (A,/?) = (90,0). The equivalent Galactic coordinates are (£,b) = (186.3725,-0.0200), so the 
corresponding velocity vector is u©, eq = (—29.814,-3.33,0.009) km s -1 . Decomposing this into source 
coordinates yields (v x ,Vy,v z )®, eq = (—10.704,19.428, —20.199) km s -1 . 

Since our model uses only a crude approximation for what in reality is a moderately noncircular orbit, 
small deviations from the model will grow with time. As such, we’d like to calibrate it using the vernal 
equinox closest in time to the observations. This occurred on 2012 May 20, which corresponds to a Modified 
Julian Date of MJD = 56006. We obtain a value of = 6.096. 
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B.3. Solar orbital motion 


The third component of the transverse motion comes from the Sun’s orbit about the Galactic center, 
relative to that of the scattering screen. From Reid et al. (2014), the distance from the Galactic center to the 
Sun is Ro = 8.34 kpc. If we denote the distance from the Sun to the scattering screen as D and the distance 
from the scattering screen to the Galactic center as R, then the law of cosines gives us an expression: 

R = sJd 2 +Rl + 2DR 0 cos(0) (B5) 

Here, 9 is the angle between £ = 180° and the direction to ESO 558-G009 (i.e., 9 = £ — 180°). 

If we define a to be the angle between the Sun and the scattering screen, as seen from the Galactic center, 
then we have a second expression for R: 

R = D cos(0 — a) + R 0 cos(a) (B6) 

Combining Equations B5 and B6 yields a numerically invertible expression for a in terms of D. Once we 
know a, we can use it to determine the component of the scattering screen’s orbital motion that lies along 
the same direction as the Sun’s orbital motion. If the orbital velocity of the scattering screen is V s , then the 
parallel component is just V s cos(a). 

The line of sight towards ESO 558-G009 is such that the scattering screen lies outside of the solar orbit. 
The rotation curve of the Milky Way is known to be very nearly flat at these outer radii (see Reid et al. 
2014), with an orbital velocity of 240 km s -1 . We can thus set V s = Vq = 240 km s -1 , and we obtain a net 
apparent motion of the scattering screen (directed along the Sun’s orbital velocity vector) of: 

Vf|=F 0 (l-cos(a)) (B7) 

The Sun’s orbital motion is directed towards the Galactic coordinates (£, b) = (90, 0), which is directed along 
the E-axis. The perpendicular component of the scattering screen’s orbital velocity (i.e., the component 
directed along the X-axis) will then just be V± = —Vq sin(a). We can now use our previously-derived unit 
vectors to transform this into the source frame. Doing so yields: 


V x = Vq 
V v = Vq 
V z = Vq 


- 0.4064sin(a) - 0.4218(1 - cos(a)) 
0.6992 sin(a) + 0.4275(1 - cos(a)) 

- 0.5882 sin(a) + 0.7996(1 - cos(a)) 


(B8) 


Combining this with Equations B3 and B4 allows us to fully characterize the transverse motion of the 
scattering screen, as seen from Earth, in terms of t (which is known for every observation) and D (which we 
would like to know). For a nearby screen, D <C Ro, and the transverse motion becomes a function of t only. 
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